function r = logmnrnd(logp)
[lsp, lcsp] = logsump(logp);
r = find(log(rand)+lsp > [-inf,lcsp], 1, 'last');

end

function [lsp, lcsp] = logsump(logp)
pivot = max(logp);
lcsp = pivot + log(cumsum(exp(logp(:)'-pivot)));
lsp = lcsp(end);

end